###8305211218_吴滔_实践5
#ANOVA方差分析####
library(reshape2);library(dplyr)
ADdata<-log2(data)#对数化
allcol<-colnames(ADdata)#提取ADdata的列名
asym<-grep(paste0("^","asym"),allcol,value = TRUE) #组内样本聚类
ad<-grep(paste0("^","ad"),allcol,value = TRUE)
ctl<-grep(paste0("^","ctl"),allcol,value = TRUE)

##单个蛋白质数据组处理
ADdata1<-ADdata[5,]#取第5行蛋白数据
group<-c(asym,ad,ctl)#创建分组
alldata<-melt(ADdata1,measure.vars = group,variable.name = "Group", value.name = "Value")#确定第一列为列名，第二列为值
alldata$Group <- gsub("\\d+", "", alldata$Group)#将Group列转换为因子类型
alldata <- alldata %>% mutate(Value = na_if(Value, -Inf))#把-Inf转变为NA，不参与计算
oneway.test(Value~Group,data = alldata)

##全部蛋白质数据处理
mylist<-list()
for (i in 1:nrow(ADdata)){
  mycount<-0#用于计数，确定有多少组观察量足够
  idata<-ADdata[i,]
  fordata<-melt(idata,measure.vars = group,variable.name = "Group", value.name = "Value")#确定第一列为列名，第二列为值
  fordata$Group <- gsub("\\d+", "", alldata$Group)#将Group列转换为因子类型
  fordata <- fordata %>% mutate(Value = na_if(Value, -Inf))#把-Inf转变为NA，不参与计算
  grouped_data<-split(fordata,fordata$Group)#组内确定样本是否足够
  for (group_name in names(grouped_data)) {
    group_data <- grouped_data[[group_name]]
    if (sum(!is.na(group_data$Value)) >= 3){mycount<-mycount+1}
  }
  if(mycount==3){mylist<-c(mylist, oneway.test(Value~Group,data = fordata)$p.value)}#每个组样本量足够才进行方差分析
  else{mylist<-c(mylist, NA)}#否则输出NA
}

#富集分析####
setwd("D:\\R_DATA")
load("volcano.RData")#加载数据源
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("clusterProfiler")
library(clusterProfiler)
BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)
BiocManager::install("enrichplot")#安装clusterProfiler、org.Hs.eg.db、enrichplot并加载
library(enrichplot)
##
proteins <- prostat[prostat$P < 0.05, "ID"]
ego <- enrichGO(gene         = proteins,
                OrgDb        = org.Hs.eg.db,  
                keyType      = "SYMBOL",      
                ont          = "BP",          # Biological Process
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05,
                readable     = TRUE)
# 绘制气泡图
bub <- enrichplot::dotplot(ego)
print(bub)
